# -*- coding: utf-8 -*-

#DESCRIPTION
#Module for computing factor of safety of slopes by various limit equilibrium methods

#Reference:
#Duncan JM and Wright SG (2005). Soil Strength and Slope Stability. John Wiley & Sons, Inc., Hoboken, New Jersey





#1. libraries and modules
#import os
#from array import *
#import csv
import HBtoMCstrength as HBMC
import math
import numpy as np
#import scipy as sp
#import scipy.optimize
#import matplotlib.pyplot as plt
#from scipy.interpolate import UnivariateSpline
#import scipy.stats.stats as st
import pandas as pd
#from datetime import datetime, date, time, timedelta

#2. math/trigo functions
tan=math.tan
sin=math.sin
atan=math.atan
rad=math.radians
pow=math.pow
exp=math.exp

#3. Limit equilibrium by Ordinary Method (Fellenius method), from Duncan and Wright, 2007
def OrdinaryMethod(s_weight,s_alpha,s_pwp,s_bl,s_SF,s_seisF,s_tcF,  
                   sign_tau,sig_t,
                   neg_sign_check, 
                   full_output_FS):

    #s_weight:          array,      slice weigth
    #s_bl:              array,      slice baselength
    #s_alpha:           array,      slice base inclination
    #s_pwp:             array,      pore water pressure at slice base
    #s_SF:              array,      driving shearing force of slice due to slice weight
    #s_tcF:             array,      force of water due to tension crack at slice
    #s_seisF:           array,      seismic force
    #sign_tau:          UnivariateSpline object
    #sig_t:             float,      minimum effective normal stress allowable
    #neg_sign_check:    boolean,    check for negative effective normal stresses if 1
    #full_output_FS:    boolean,    return arrays of parameters for all slices if 1

    #1. computing effective normal force, N1, and stress, n1
    N1=s_weight*np.cos(s_alpha) - s_pwp*s_bl*np.cos(s_alpha)**2     #eq. 6.57
    n1=N1/s_bl      #eq.6.57

    #2. checking for negative normal stresses
    if neg_sign_check==1:
        if len(np.where(n1<sig_t)[0])>0:
            return -120,[],[],[]    #returns error code -120 for any occurrence of negative normal stresses less than the material tensile strength
               
    #3. computing phi and coh from effective normal stresses
    phi, coh=HBMC.phi_coh_from_HoekBrownFE_spline(n1, sign_tau, sig_t)                          

    #4. computing shear strength, check for negative shear strength
    Rf=np.where(np.isfinite(coh),(coh+n1*np.tan(phi))*(s_bl),np.zeros(len(coh)))   #eq. 6.46

    #5. check for negative shear strength
    Rf=np.where(Rf>0,Rf,np.zeros(len(Rf)))

    #6. computing potential force of water in tension crack 
    if len(np.where(n1<=0)[0])>0:
        tcF=s_tcF[np.where(n1<=0)[0]]
    else:
        tcF=0

    #7. computing factor of safety
    FS=np.sum(Rf)/(np.sum(s_SF)+ np.sum(s_seisF))#+ tcF)       #6.48, 6.73

    #8. preparing data for output 
    if full_output_FS==1:    
        Fel=pd.DataFrame(data=None,index=range(1,len(n1)+1))
        Fel['N1']=N1
        Fel['n1']=n1
        #Fel['n1pos']=n1pos
        Fel['phi']=phi
        Fel['coh']=coh
        Fel['RF']=Rf
        Fel['SF']=s_SF
        Fel['tcF']=s_tcF
        Fel['seisF']=s_seisF
    else:
        Fel=[]
    
    return FS,phi,coh,Fel

#4. Limit equilibrium by Simplified Bishop method, from Duncan and Wright, 2007
def SimplifiedBishop(FS,
                     phi,coh,s_weight,s_bl,s_alpha,s_pwp,s_SF,s_tcF,s_seisF,
                     sign_tau, sig_t,
                     malpha_check, neg_sign_check,converge_check,
                     full_output_FS):

    #FS:                float,      initial estimate of factor of safety, usually the result of analysis by Ordinary Method
    #phi, coh:          array,      phi and cohesion value of slices from previous analysis by Ordinary method
    #s_weight:          array,      slice weigth
    #s_bl:              array,      slice baselength
    #s_alpha:           array,      slice base inclination
    #s_pwp:             array,      pore water pressure at slice base
    #s_SF:              array,      driving shearing force of slice due to slice weight
    #s_tcF:             array,      force of water due to tension crack at slice
    #s_seisF:           array,      seismic force
    #sign_tau:          UnivariateSpline object
    #sig_t:             float,      minimum effective normal stress allowable
    #malpha_check:      boolean,    check m-alpha if 1
    #neg_sign_check:    boolean,    check for negative effective normal stresses if 1
    #converge_check:    boolean,    check solution convergence if 1
    #full_output_FS:    boolean,    return arrays of parameters for all slices if 1

    #1. initializong values for iterative solving
    FS_diff=100
    max_iter=50
    itr=0

    #2. iterative solving
    while FS_diff>0.005:

        #3. computing normal force, N, and effective normal stress, n1
        N=((s_weight - (1/FS)*(coh*s_bl-s_pwp*s_bl*np.tan(phi))*np.sin(s_alpha))/
        (np.cos(s_alpha)+np.sin(s_alpha)*np.tan(phi)/FS))   #eq. 6.77

        n1=N/(s_bl)-s_pwp    #eq. 6.65

        #4. computing phi and cohesion from effective normal stresses
        phi, coh=HBMC.phi_coh_from_HoekBrownFE_spline(n1, sign_tau, sig_t)                         

        #5. computing malpha
        malpha=np.cos(s_alpha)*(1 + np.tan(s_alpha)*np.tan(phi)/FS)    #eq. 6.77

        #6. computing shear strength, check for negative shear strength
        Rf=np.where(np.isfinite(coh),(coh+n1*np.tan(phi))*(s_bl),np.zeros(len(coh)))   #eq. 6.46

        #7. check for negative shear strength
        Rf=np.where(Rf>0,Rf,np.zeros(len(Rf)))

        
        #8. computing potential force of water in tension crack 
        if len(np.where(n1<=0)[0])>0:
            tcF=s_tcF[np.where(n1<=0)[0]]
        else:
            tcF=0

        #9. computing factor of safety
        FS_new=np.sum(Rf)/(np.sum(s_SF)+ np.sum(s_seisF))#+ tcF)       #6.78,
    
        #10. processing parameters for iteration
        FS_diff=abs(FS-FS_new)
        FS=FS_new
        itr=itr+1
        if itr>max_iter:break

    #11. preparing data for output 
    if full_output_FS==1:
        Bis=pd.DataFrame(data=None,index=range(1,len(n1)+1))
        Bis['N']=N
        Bis['n1']=n1
        Bis['coh']=coh
        Bis['phi']=np.degrees(phi)
        Bis['malpha']=malpha
        Bis['RF']=Rf
        Bis['SF']=s_SF
        Bis['tcF']=s_tcF
        Bis['seisF']=s_seisF
    else:
        Bis=[]

    #12. check for iteration convergence
    if converge_check==1:
        if itr>max_iter:
            return -111,Bis #returns error code -111 for non-convergent solution

    #13. check for m-alpha<0.2
    if malpha_check==1:
        if len(np.where(malpha<=0.2)[0])>0:  
            return -112,Bis #returns error code -112 for very steep slip surfaces

    #14. check for negative effective normal stress
    if neg_sign_check==1:
        if len(np.where(n1<sig_t)[0])>0:
            return -120,Bis #returns error code -120 for any occurrence of negative normal stresses less than the material tensile strength

    return FS,Bis



        

        

            
            
            

                

    











    
